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ABSTRACT 
Two  stationary  first-order  autoregressive  processes  with  Beta  marginal 
distributions  are  presented.  They  are  both  linear,  additive  processes  but  the 
coefficients  are  Beta  random  variables.  Their  autocorrelation  functions  are 
investigated:  one  is  positive  and  the  other  alternates  in  sign.  The  useful- 
1  ness  of  the  models  in  simulation  is  discussed.  The  Bivariate  Beta  distribu- 
tions of  two  consecutive  observations  are  considered  in  some  detail.  Several 
examples  are  given,  including  a  Bivariate  Uniform  process  which  is  also 
examined  in  detail.  The  relationship  of  these  Bivariate  Beta  distributions 
to  the  Dirichelet  distribution  is  discussed. 


1 

1.  INTRODUCTION 

The  Beta  distribution  is  the  most  versatile  and  useful  distribution 
available  on  a  bounded  interval.  Despite  this  there  are  very  few  practical 
models  for  describing  correlation  between  pairs  of  Beta  random  variables  or 
serial  correlation  in  sequences  of  them.  This  is  unfortunate  since  there  is  a 
natural  interest  in  modelling  sequences  of  dependent  Beta  variates.  These 
arise  in  a  variety  of  ways  but  are  often  associated  with  the  stochastic 
behaviour  of  a  proportion  or  probability  over  time.  One  example  is  the  study 
of  market  share,  e.g.  the  proportion  of  the  market  held  by  a  particular 
product.  (See,  for  example,  Wichern  and  Jones  (1977)). 

Some  efforts  have  been  made  to  model  such  behaviour.  Azzalini  (1982) 
developed  a  simple  Markov  model  for  use  in  a  quality  control  context  where  the 
proportions  defective  in  adjacent  batches  sampled  may  be  expected  to  be  depen- 
dent. He  used  the  Product  Autoregressive  process  discussed  by  McKenzie  (1983). 
This  model  has  limited  usefulness  for  the  Beta  distribution,  however,  since  it 
requires  that  of  the  two  parameters  of  that  distribution  one  must  be  integral 
and  the  other  rational.  Another  approach  has  been  suggested  by  Souza  and 
Harrison  (1981).  They  attribute  to  the  proportion  a  Beta  distribution  which  is 
revised  with  each  observation  by  a  procedure  based  on  Bayesian  and  Information 
Theoretic  concepts.  This  approach  appears  to  offer  some  promise  for  forecast- 
ing but  is  fairly  complex  in  general.  A  more  traditional  time-series  approach 
to  such  a  data  set  would  be  to  treat  the  proportions  X   not  as  Beta's  but  as 
Logistic-normal  variates.  Thus,  we  would  transform  X   to  log-odds,  i.e. 
Y  =  £n[X  /(1-X  )]  ,  and  attempt  to  model  {Y  }  as  an  autoregressive  moving- 
average  process.  This  approach  is  suggested  by  the  work  of  Aitchison  and  Shen 
(1980)  and  Aitchison  (1982),  but  does  not  appear  to  have  been  investigated  as 
a  time-series  procedure  yet.  It  discards  the  Beta  marginal  distribution  which 
is  our  main  interest  here.  Further,  it  appears  to  be  extremely  difficult  to 


reverse  the  procedure  and  generate  a  sequence  of  proportions  with  a  specific 
correlation  structure. 

The  purpose  of  this  work  is  to  present  a  discrete  time  Markov  process  with 
a  Beta  marginal  distribution.  It  is  constructed  in  the  spirit  of  the  recent 
work  on  modelling  of  non-Gaussian  time  series  illustrated,  for  example,  by 
Lawrance  and  Lewis  (1980,  1982)  and  Jacobs  and  Lewis  (1983).  Since  simulation 
is  also  a  major  motivation  for  such  work,  we  seek  models  which  are  simple  and 
flexible  and  whose  parameters  are  few  and  physically  meaningful.  The  aspect  of 
simulation  is  important  here  for,  as  Schmeiser  and  Lai  (1980)  noted  in  a  recent 
survey  (1980),  there  are  few  practical  ways  of  generating  dependent  Beta  random 
variables.  This  is  because  the  usual  multivariate  Beta  distributions  are 
closely  related  to  the  Dirichelet  distribution  and  so  constrain  their  vector 
variates  in  a  way  which  is  undesirable  in  general. 

We  present  here  two  simple  discrete-time  stochastic  processes  whose  mar- 
ginal distributions  are  Beta  random  variables.  The  processes  are  linear, 
additive  autoregressions  with  random  coefficients.  The  coefficients  themselves 
are  also  Beta  variates.  There  is  a  single  free  parameter  which  corresponds  in 
a  simple  way  to  the  correlation  in  each  model  and  a  wide  range  of  correlation 
is  possible.  Because  of  their  simplicity,  the  models  provide  a  powerful  way  of 
generating  sequences  of  dependent  Beta  variates  using  only  independent  Betas. 
As  noted  above,  there  is  a  scarcity  of  practical  multivariate  Beta  distribu- 
tions. Thus,  the  bivariate  distributions  associated  with  the  processes  are 
discussed  in  some  detail.  They  exhibit  a  number  of  interesting  features.  We 
also  examine  in  detail  the  particular  case  of  the  bivariate  Uniform  distribu- 
tion. It  plays  an  important  role  in  the  simulation  of  dependent  pairs  of 
random  variables . 


Before  describing  the  models  we  may  not  that  there  is  a  well-known 
continuous-time  Markov  process  with  Beta  marginals.  It  arises  in  genetics  and 
is  one  of  the  forms  of  the  Wright-Fisher  gene  frequency  models  described  in 
detail  in  the  books  by  Karlin  and  Taylor  (1975,  1981).  It  is  a  diffusion 
process  and,  in  different  forms,  has  found  applications  in  sociology,  psychol- 
ogy and  marketing.  It  is  also  derived  by  Massey  et  al .  (1970)  as  a  stochastic 
response  model.  They  obtain  it  as  the  limiting  form  of  the  "contagious 
binomial"  distribution  developed  by  Coleman  (1964)  to  model  voting  behaviour. 
Of  course,  such  processes  are  in  continuous  time  and  it  is  by  no  means  clear 
how  they  can  be  restructured  in  discrete  time.  Nor  is  it  clear  whether  the 
processes  presented  here  represent  some  discrete  time  formulations  of  the 
diffusion  processes. 

2.  THE  MODELS 

2.1.  A  random  variable  (r.v.)  X  is  said  to  have  a  Beta  distribution  with 

parameters  (a, 3)  if  it  has  probability  density  function  (p.d.f.) 

f(*> •    B(i;x6r  ■  ° < x < i s «.»><> • 

For  convenience  in  what  follows,  we  shall  write  such  a  random  variable  as 
Be(a,3)  .  We  note  for  later  use  that  for  such  X  , 

E(X)  =  a/(a+B),  Var(X)  =  a6/[( a+$)2( a+8+1) ]  and  the  third  moment  about  the 
mean  is  given  by  nu  =  2a6(  B-a)/[(a+B)  ( a+B+1)  ( a+0+2)  ]  . 
The  models  presented  here  use  the  following  results: 

1  -  Be(a,3)  =  Be(e,a)  (1) 

Be(a,3)  •  Be(a+6,y)  =  Be(a,B+y)  (2) 


The  first  of  these  two  results  is  well  known  and  easily  demonstrated.  The 
second  result  states  that  the  product  of  two  independent  Beta  r.v.s.  with 
parameters  as  specified  is  itself  a  Beta  r.v.  The  result  may  be  verified  by 

c 

considering  the  Mel  Tin  Transform  (Widder,  1946),  i.e.  E(X  )  ,  of  the  product 
on  the  left-hand  side  of  (2).  It  is 

r/yS^  •  B(<*+S>g)  .  B(g+e+S,y)  _  TJ  g+S )  T{  ci+3+y)  _  B(a+S,3+y) 
LlA  ;  ~  B(a,3)     B(a+3,Y)    r(  a)  r(  a+g+y+s)   B(a,3+y) 

which  is  the  transform  of  the  right-hand  side  of  (2).  We  shall  refer  to  the 
application  of  (2)  to  change  one  Beta  r.v.  into  another  as  the  Beta-Beta 
transformation. 

Two  distinct  models  are  presented  here.  One  is  for  positively  correlated 
pairs  of  Beta  r.v.s.  and  is  denoted  by  PBAR  and  the  other  is  for  negatively 
correlated  pai^s  of  Beta's  and  is  denoted  by  NBAR.  Both  models  are  linear  and 
additive  and  have  random  coefficients. 

2.2.  The  PBAR  model  is  given  by 

Xt  =  1  -  Ut(l  -  »tXt_x)  (3) 

where  {U.}  and  {W  }  are  independent  sequences  of  independent  identically 
distributed  (i.i.d.)  r.v.s.,  independent  of  previous  X's  ,  and  U   is 
Be(3,a-p)  and  W   is  Be(p.a-p),  (0  <  p  <  a)  .  Now  if  X  ,   is  Be(a,3)  then 
W .X  .   is  Be(p,a+3-p)  from  (2)  and  1  -  W  X  ,   is  Be(a+3-p,p)  from  (1). 
Further  use  of  (2)  shows  that  U  (1-W  X  ,)  is  Be(3,a)  and  so  X   given  by  (3) 
is  Be(a,3)  by  (1).  Thus,  equation  (3)  defines  a  stationary  process  {X  }  with 
a  Beta(a,3)  marginal  distribution. 

Further,  there  is  a  single  free  parameter  in  this  scheme,  viz.  p  .  As  we 
shall  see,  the  value  of  p  determines  the  correlation  structure  of  the 


process  {Xt}  .  From  the  structure  of  (3),  the  process  is  a  first-order  linear 
autoregression  with  random  coefficients.  Direct  calculation  yields  the  auto- 
correlation  function  of  the  process  as  px(k)  =  p  ,  k  =  0,1,...,  where 

p  =  E(U)E(W)  =  pe/a(a+6-p)  .  (4) 

Now  p  as  defined  by  (4)  is  a  monotonic  increasing  function  of  p  for  fixed 
(a, 6)  .  Further,  since  e  >  0  ,  0  <  p  <  a  ,  we  may  deduce  that  0  <  p  <  1  . 
Thus,  the  entire  range  of  positive  correlation  is  possible  for  any  values  of  a 
and  3  ,  i.e.  any  Beta  marginal  distribution. 

There  are  two  limiting  cases  for  the  parameter  p  we  may  consider.  When 
p  is  zero  W  =  0  with  probability  one  and  U,   is  Be(0,a)  .  Thus,  X   is 
independent  of  X  .   and  p  =  0  .  When  p  =  a  both  U.   and  W   are  unity 
with  probability  one.  Thus,  X  =  X.  .  and  p  =  1  .  The  process  is  not  ergodic 
in  this  case. 

It  is  important  to  notice  that  with  this  model  (3)  p  =  0  implies 
independence. 

2.3.  The  model  for  the  NBAR  process  is  given  by 

xt =  V1  -  wtxt-i»  •  (5) 

As  before,  {V  }  and  {W  }  are  independent  sequences  of  i.i.d.,  r.v.'s 
independent  of  X  ,   and  V   is  Be(a,3-p)  and  W   is  Be(p,a-p)  . 
Again,  it  is  easily  verified  using  (1)  and  (2)  that  equation  (5)  will  generate 
a  stationary  process  whose  marginal  distribution  is  Beta(a,e)  .  The  NBAR 
process  is  also  a  linear  autoregression  with  random  coefficients.  It  too  has 
autocorrelation  function  of  the  form  p  ,  k  =  0,1,...,  but  now 

p  =  -  p/(a+6-p)  •  (6) 


Notice  that  the  specification  of  the  distributions  of  V   and  W   requires 
that  0  <  p  <  a  and  0  <  p  <  0  .  From  (6),  p  is  a  monotonic  decreasing 
function  of  p  for  fixed  (a, 3),  and  so  we  may  deduce  that  for  the  NBAR 
process(5)  -max(a/3,3/a)  <  p  <  0  . 

The  upper  extreme  of  zero  is  again  attainable  when  p  =  0  and  X  =  V  , 
which  is  independent  of  X.  ,  .  As  with  the  PBAR  process,  it  is  important  to 
note  that  p  =  0  implies  idependence.  The  lower  limit  -max(a/3,3/cc)  is  also 
attainable.  If  3  <  a  it  is  attained  when  p  =  3  and  so  V.  =  1  .  If 
3  >  a  it  occurs  where  p  =  a  which  corresponds  to  W  =  1  .  When  a  =  3 
the  lower  limit  is  -1  which  corresponds  to  the  usual  antithetic  relationship, 
Xt  =  1  -  Xt_1  ,  given  by  Vt  =  Wt  =  1  . 

2.4.  The  models  PBAR  and  NBAR  given  by  (3)  and  (5)  yield  random  coefficient 
autoregressions  of  order  one.  Further,  the  first-order  correlation  p  satisfies 

-  max(|  ,  £)  <  P  <  1  .  (7) 

This  is  not  the  greatest  possible  range  for  general  (a, 3)  .  For  example,  if 
a  =  2 ,  3  =  1  we  find  that  -0.5  <  p  <  1  for  these  models.  On  the  other  hand, 
it  may  be  deduced  from  Moran  (1967)  that  the  correlation  between  two  Be(2,l) 
r.v.s.  is  bounded  below  by  (9tt-32)/4,  i.e.  approximately  -0.9314.  How  far  the 
lower  bound  given  by  (7)  is  from  the  minimum  correlation  possible  is  not  known  for 
general  (a, 3)  . 

We  may  note,  however,  that  in  the  symmetric  case,  i.e.  a  =  3  ,  the  range 
is  (-1,1)  as  would  be  hoped.  As  before,  the  two  extremes  may  be  attained: 
the  upper  limit  of  1  from  X.  =  X.  ,   and  the  lower  limit  of  -1  from  the 
usual  antithetic  relationship  X.  =  1  -  X.  ,  .  This  latter  is  obtained  from 
NBAR  with  p  =  a  =  3  so  that  V  =  Z  =  1  .  In  particular,  it  is  now 


possible  to  generate  r.v.s.  uniform  on  (0,1)  with  any  first-order  autoregressive 
correlation,  i.e.  any  pe[-l,l]  .  We  return  to  this  point  later. 

3.  BIVARIATE  DISTRIBUTIONS  OF  (X  ,  X   ) 

3.1.  For  convenience,  we  rewrite  the  PBAR  model  (3)  as  Y  =  1  -  U(l-WX)  . 
Initially,  consider  the  joint  p.d.f.  of  (Y,X)  conditional  on  U  .  It  can  be 
written  in  the  form 

fY,x|u"'xlu>  'h-  ¥*>  •  V^'  •     i-y<"<&- 

The  joint  p.d.f.  of     Y,  X     and     U     may  now  be  derived  and  from  it  the  p.d.f. 
of  (Y,X)    is  obtained  in  the  form 

, .      v  3-1  m(y  ,x)        ..  -  . 

ft  x(y,x)   =  ii=il       /  sp-1[x(l-y)-(l-x)s]a-P-1(l-y+s)e-a(y-s)a^-1ds  (8) 

where  m(y,x)  =  min(y  ,x( 1-y)/ (1-x) )  and  C+  =  B(a,e)B(p,a-p)B( 6,a-p)  . 

A  change  of  variable  t  =  (l-x)s/(l-y)  in  the  integral  in  (8)  yields  the 
result  fv  Y(y>x)  =  ^v  v(x»y)  •  Thus,  the  bivariate  p.d.f.,  which  is  defined 

Y  ,  A  I  ,  A 

on  the  unit  square,  is  symmetric  about  y  =  x  .  Note  also  that  m(y,x)  =  y 

if  y  <  x  .  Thus,  the  two  forms  of  the  integral  exist  on  either  side  of 

y  =  x  ,  and  we  can  define 

g(y,x)     y  <  x 

fy>x(y,x)  =  (9) 

g(x,y)     y  >  x 

where  g(y,x)  is  given  by  the  right  hand  side  of  (8)  with  m(y,x)  replaced 
by  y  . 


This  symmetry  is  also  important  from  the  viewpoint  of  modelling  or  iden- 
tification of  the  PBAR  process.  The  importance  arises  from  the  idea  of  time- 
reversibility  of  a  stationary  process.  The  concept  is  discussed  in  detail  by 
Weiss  (1975).  A  discrete-time  stationary  process  {X.}  is  time-reversible  if 
the  joint  distributions  of  {X-  ,X?,. ..  ,X  }  and  {X  ,X  l,..,,X-}  are  identical 
for  every  t  .  In  the  case  of  a  first-order  autoregression,  as  here,  the 
joint  p.d.f.  can  be  expressed  as 

t-1  t-1 

f(x.  ,x?,...,x  )  =  n  f(x.,x.  .)/  n  f(x.) 
i  c  t    i=1    111   i=1    i 

Such  a  process  is  time  reversible  then  whenever  the  bivariate  distribution  is 
symmetric.  Thus,  the  PBAR  process  is  time-reversible. 
Indeed,  the  time  reversed  process  is  yiven  by 

xt  - i  -  u;(i-w;w 

where  {U'  }  and  {WM  have  exactly  the  same  properties  as  {U.}  and  {W.} 
respectively  defined  by  (3). 

3.2.  For  the  NBAR  process,  writing  Y  =  V(l-WX)  and  proceeding  as  for  the  PBAR 

process  yields 

n   n  3-1  m(l-y,x)   ,  i      o   i 

fY,X(y'x)  =  iTI /       s   [xy  -  (l-xJsl^P^d-y-sJ^P^ds       (10) 

where  C_  =  B(a,S)B(p,o-p)B(a,8-p)  . 

Again,  the  change  of  variable  t  =  (l-x)s/(l-y)  shows  that  the  bivariate 
p.d.f.  is  symmetric  about  y  =  x  and  the  NBAR  process  is  time-reversible. 
The  structure  of  the  density,  however,  is  a  little  more  complex  than  for  the 
PBAR  process.  The  upper  limit  of  the  integral  in  (10)  is  m(l-y,  x)  =  1  -  y 


if  x  +  y  >  1  .  Thus,  the  form  of  f~  depends  upon  which  side  of  the  line 

x  +  y  =  1  we  are  on. 

Define 

9i(y,x)      x  +  y  <  1  , 

fY  x(y»x)  =  (ID 

92(y,x)      x  +  y  >  1  . 

From  symmetry  about  y  =  x  ,  we  know  that  g,(y,x)  =  g,(x,y)  and 
92(y»x)  =  9p(x>.y)  *  °^  more  immedlate  interest  is  the  relationship  between 
g   and  g~  .  A  further  change  of  variable  in  the  integral  in  (10)  yields  the 
following  result.  Using  an  obvious  notation  to  denote  dependence  upon  the 
parameters  of  the  marginal  distribution: 

g^l-y,  1-x;  a,  3)  =  g2(y,  x;  3,  a)  (12) 

Using  (12)  to  evaluate  f"  in  both  the  triangular  regions  induced  by  the  line 
x  +  y  =  1  yields 

f~  x(l-y,l-x;  ot,e)  =  fy  x(y,x;  6, a)  (13) 

This  result  (13)  is  an  obvious  two-dimensional  analogue  of  the  relationship 

specified  by  equation  (1),  viz.  fw(l-x;  a, 3)  =  'M*'  e'a^  * 

Further,  by  the  symmetry  about  y  =  x  ,  equation  (12)  yields 

g1(l-x,  1-y;  a,e)  =  g2(y,x;  B,a)  (14) 

which  specifies  the  nature  of  the  relationship  between  g.  and  g„  across 
the  line  x  +  y  =  1  .  In  particular,  note  that  when  a  =  3  the  bivariate 
p.d.f.  is  symmetric  about  both  y  =  x  and  x  +  y  =  1  . 
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4.  MOMENTS 

Although  it  is  difficult  to  obtain  the  p.d.f.s  in  an  explicit  form,  much 
information  can  be  derived  directly  from  the  structural  relationships,  (3)  and 
(5).  This  is  also  true  of  the  conditional  distributions.  In  both  cases,  the 
conditional  expectation  is  linear.  For  the  PBAR  process,  using  (3)  yields 

E(Xt+1|Xt  =  x)  =  [(a+B-p-ae)  +  P0x]/a(  a+3-p)  ,  (15) 

and  for  the  NBAR  process,  (5)  yields 

E(Xt+1|Xt  =  x)  =  (a-px)/(o+B-p)  .  (15a) 

Further,  the  time-reversibility  of  the  processes  ensures  that  the  inverse 
regressions  are  identical  to  (15)  and  (15a),  i.e.  E( X. | X.  ,  =  x)  =  E(X.+, |X.  =  x) 

In  both  cases  the  conditional  variances  are  quadratic  and  of  the  form 

2     2     2  2  2  2         2 

var(XJ.,1|X.  =  x)  =  a  (1-mx)  +  (a  +m  )a,,x   where  a, ,  =  var(W)  and  m  and 
t+1 '  t  w  W  W 

2 
a       are  the  mean  and  variance  of  U  for  PBAR  and  of  V  for  NBAR. 

Higher  order  moments  are  important  in  the  identification  of  non-standard 

time-series  models  and  we  note  here  two  of  particular  interest.  The  first  is 

Cp,(k)  =  Cov(X.,  X.  ,)a.     For  time  reversible  processes  C?,(k)  =  C«,(-k)  . 

i/ 
For  both  the  PBAR  and  NBAR  processes,  C?,(k)  =  m~p  ,  k  =  0,1,2,...,  where 

m~  is  the  third  moment  of  X   about  its  mean,  y  say,  and  these  are  given 

in  Section  2.1 . 

Another  moment  of  particular  value  in  residual  analysis  is  based  on  the 

residuals  from  minimum  mean  square  error  prediction,  i.e. 

Rt  =  (X.  -  u)  -  p(X.  ,  -  u)  •  Such  residuals  are  uncorrelated  for  PBAR  and 

2 
NBAR  processes.  Of  more  interest  is  the  behaviour  of  Cov(R?,  R  ,  )  which 

is  useful  in  distinguishing  between  constant  and  random  coefficient  models. 

For  both  the  PBAR  and  NBAR  processes 
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2  k 
m3(l-p)(l-p  )p  , 


Cov(fV  Rt-k)  =  ° 


m3(l+2p)(l-P)\ 


k  =  1,2,.. 
k  =  -1,-2, 
k  =  0  . 


Further  details  and  discussion  of  the  usefulness  of  these  moments  may  be  found 
in  Lewis  and  Lawrance  (1983b). 


5.  EXAMPLES 

To  illustrate  the  nature  of  the  bivariate  p.d.f.  the  functions  g  from 
(9)  and  g,  and  g?  from  (11)  are  evaluated  explicitly  for  a  few  specific 
cases  below. 
PBAR 
PI:  a  =  3  =  2;  p  =  1  i .e.  p  =  1/3, 

g(y,x)  =  12y(l-x) 
P2:  a  =  3,  6  =  2;  p  =  1  i.e.  p  =  1/6, 

g(y,x)  =  72(l-x)[(l-x)y(2-y)-2(l-yUn(l-y)-2y(2-x-y)] 
P3:  a  =  4,  3  =  2;  p  =  3  i.e.  p  =  1/2, 


g(y,x)  =  120(l-x)[l+2(l-y)*n(l-y)  -  {l-yfl 


NBAR 
Nl: 


a  - 


■  2;  p  =  1 

g1(y,x 


g2(y»x 

N2:   a=3,B=2;p= 

92(y,x 

N3:  a  =  5,  B  =  13/3;  p  =  4  i.e.  p  :    3/4 

g2(y,* 


g^y.x 


.e.  p  =  -  1/3, 
■  12xy 
=  12(l-x)(l-y) 

1  i.e.  p  =  -  1/4 

,,2  2 
=  36x  y 

=  36(l-x)(l-y)(x+y  +  xy-1) 


=  A(l-x)10/3(l-y)10/3 


=  g2(y,x)  -  B(l-x-y)1/3{3(l-y)3(l-x)3  -  i(  l-y)2( l-x)2( 1-x-y) 

+  y(l-y)(l-x)(l-x-y)2  -  ^(l-x-y)3} 
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where  A  =  3532100/2187  =  1615.04,  and  B  =  140A/283  =  798.96. 
Plots  of  the  contours  and  the  surface  of  PI  are  displayed  in  Figures  1(a)  and 
(b).  If  they  are  rotated  through  90°  about  the  point  (0.5,  0.5)  the  corre- 
sponding plots  are  obtained  for  Nl .  Contour  plots  are  displayed  for  P2  and  N2 
in  Figures  2  and  3. 

Figures  1 ,2,3 

The  most  remarkable  feature  of  these  distributions  is  the  appearance  of  a 
ridge  in  the  surface  of  PI  (and  so  Nl)  but  not  in  P2  or  N2.  The  ridge  is  due 
to  the  fact  that  both  densities  have  two  forms  as  described  by  equations  (9) 
and  (11).  Thus  the  ridge,  if  it  occurs,  corresponds  to  the  line  y  =  x  for 
PBAR  p.d.f.s  and  x  +  y  =  1  for  NBAR's.  From  the  definition,  the  p.d.f.s  will 
be  continuous  but  tnei r  derivatives  need  not  be  so.  Any  points  of  disconti- 
nuity occur  on  these  two  lines.  To  determine  the  general  conditions  for 
occurrence  of  a  ridge  we  can  examine  the  behaviour  of  the  derivatives  of  the 
p.d.f.s  on  both  sides  of  the  two  lines.  This  procedure  yields  the  following 
results.  For  the  PBAR  density  (8)  there  is  no  ridge  provided  p  <  (a-1)  . 
Using  (4),  we  may  deduce  that  no  ridge  occurs  provided  p  satisfies 

0  <_  p  <  0(o-l)/a(0+l)  .  (16) 

For  the  NBAR  density  (10)  there  is  no  ridge  provided  p  <  (a+3-2)/2  i.e. 
provided 

0  >  p  >  [2  -  (a+B)]/(2+cd-8)  .  (17) 

These  conditions  (16)  and  (17)  are  violated  by  P3  and  N3  respectively  and  so 
both  densities  have  a  ridge.  These  are  illustrated  in  the  surface  plots  of 
P3  and  N3  displayed  in  Figures  4  and  5,  respectively. 

Figures  4  ,5 
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6.  THE  UNIFORM  PROCESS 

A  bivariate  distribution  of  particular  interest  is  that  with  Uniform 
marginals.  Apart  from  the  natural  interest  in  modelling  data  from  such  a  dis- 
tribution it  has  an  important  application  in  simulation.  By  using  the  inverse 
distribution  function  transformation  of  the  Uniforms  bivariate  distributions 
with  any  other  marginals  can  be  obtained.  This  is  an  important  approach  to 
the  generation  of  pairs  of  dependent  random  variables. 

Some  recent  development  of  Bivariate  Uniform  distributions  appears  in  the 
papers  by  Barnett  (1980)  and  Lewis  and  Lawrance  (1983a).  The  former  exhibits 
several  different  Bivariate  Uniform  distributions  but  they  are  generally 
complex,  have  limited  correlation,  and  require  distribution  function  trans- 
formations to  obtain  the  Uniforms.  The  latter  work  gives  several  procedures 
for  generating  a  pair  of  dependent  Uniforms  from  a  random  coefficient  regres- 
sion on  independent  Uniforms.  The  procedures  are  simple  and  the  entire  range 
of  correlation  can  be  attained.  However,  to  achieve  this  breadth  the  correla- 
tions are  usually  complex  functions  of  two  parameters. 

The  Bivariate  Uniform  with  p.d.f.s  given  by  (8)  and  (10)  is  particularly 
useful  in  this  kind  of  application.  The  entire  range  of  correlation  [-1,  1]  is 
available  to  the  Uniforms  and  so  the  entire  possible  range  will  be  available 
to  the  transformed  variates.  Further,  the  correlation  is  a  simple  function  of 
a  single  parameter,  i.e.  p  =  ±  p/(2-p)  ,  (0  <_  p  <_  1)  and  so  any  desired  cor- 
relation is  easily  achieved.  Finally,  the  generation  of  the  desired  pair  is 
straightforward  involving  only  the  additional  generation  of  two  independent 
Beta  r.v.s. 

Since  the  Bivariate  Uniform  corresponds  to  (8)  and  (10)  with  a  =  e  =  1 
considerable  simplification  is  possible.  In  this  case, 

fv  y(y»x)  =  fv  v(l-y»x)  and  so  on^  f   is  considered  here.  Making  a  suit- 
able  change  of  variable  in  the  integral  yields  a  more  useful  expression  for 
the  density  given  by  (9),  viz. 
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B(p,l-p)xp(l-y)p  0  xu  yj 

The  density  has  a  line  singularity  on  y  =  x  as  can  be  seen  from  the  contour 
and  surface  plots  shown  in  Figure  6  for  the  case  p  =  1/2  . 

Figure  6 

7.  RELATIONSHIP  WITH  THE  DIRICHELET  DISTRIBUTION 

The  best-known  "Bivariate  Beta"  distribution  is  the  Dirichelet  distrib- 
ution. It  is  described  in  detail  by  Johnson  and  Kotz  (1972).  It  is  defined 
not  on  the  Unit  square  as  are  the  densities  given  by  (8)  and  (10)  but  on  the 
Unit  simplex,  i.e.   {(x,y):  0<x,y<l,x+y<l}.  As  such  it  plays  a 
natural  role  as  the  joint  distribution  of  two  proportions  from  a  single  popu- 
lation. It  seems  to  be  generally  regarded  as  a  Bivariate  Beta  because  both 
the  marginal  and  conditional  distributions  are  Beta.  The  joint  p.d.f.  in  the 
case  of  identical  Be(a,B)  marginals  is 

f(u  ^  _    ,r(g+B)      .   a-1  o-ln    ^8-a-l 

f(y>x)  "  r(o)r(a)r(s-a)   x  y   (1  y~x) 

and  the  correlation  between  X  and  Y  is  p  =  -  a/B  (a  <  3)  .  Note  that  p 
here  depends  explicitly  upon  the  parameters  of  marginal  distribution  and  so  is 
fixed  for  any  particular  marginal  distribution.  Further,  -a/B  is  the  minimum 
correlation  attainable  from  the  NBAR  process.  As  noted,  it  is  attained  when 
p  =  a  and  Y  =  V(l-X)  where  V  is  Be(a,B-a)  and  independent  of  X  . 
Thus,  Y/(l-X)   is  independent  of  X  ,  and  symmetry  ensures  that  X/(l-Y)   is 
independent  of  Y  .  This  characterizes  (X,Y)  as  having  a  Dirichelet  distrib- 
ution by  a  result  of  Darroch  and  Ratcliff  (1971).  Thus,  we  may  view  the 
Dirichelet  distribution  as  a  limiting  form  of  the  NBAR  process  distribution, 
as  p  ^  min(a,B)  . 
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8.  EXTENSIONS 

It  is  possible  to  extend  the  time-series  models  (3)  and  (5)  to  higher 
orders  of  dependence.  However,  two  simple  and  more  immediate  extensions  lie 
closer  to  the  area  of  simulation  and  are  noted  here.  The  fact  that  the  PBAR 
and  NBAR  models  yield  simple  but  powerful  methods  of  generating  sequences  of 
correlated  Beta  r.v.s.  has  been  emphasized.  Such  sequences  are  stationary  so 
that  each  Beta  r.v.  has  the  same  distribution.  An  obvious  extension  is  to  the 
generation  of  pairs  of  dependent  Beta  variates  with  different  distribution. 
This  may  be  achieved  in  a  variety  of  ways  using  (1)  and  (2)  and  (3)  or  (5)  if 
we  wish. 

A  second  simple  generalization  is  to  bounded  intervals  other  than  (U,l). 
Since  an  alternative  sample  space  is  achieved  by  a  linear  transformation  the 
procedure  is  straightforward  and  all  correlations  are  unaffected.  Thus,  sup- 
pose we  wish  to  develop  a  PBAR  process  for  Be(a,B)  r.v.s.  defined  on  (a,b) 
rather  than  (0,1).  By  considering  the  X. ' s  transformed  to  Y.'s  on  (0,1) 
and  {Y  }  satisfying  (3)  we  find  that  the  PBAR  for  {X  }  is  given  by 

X  =  b  -  U  [b-a-Wt(Xt-a)]  . 
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Figure  1(a).  Bivariate  Beta  P1,  p=0.33,  contours  .5(.5)2.5 
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Figure  1(b).  Bivariate  Beta  PI.  p=0.33.  Surface  Plot 
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Figure  3.  Bivariate  Beta  N2,  p=-0.25,  contours  .5(.5)3 
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Figure  4.  Bivariate  Beta  P3,  p  =  0.5 
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Figure  5.  Bivariate  Beta  N3,  p  =  -0.75 
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